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We study the ground state phases of the 5=1/2 Heisenberg quantum antifer- 
romagnet on the spatially anisotropic triangular lattice and on the square lattice 
with up to next-next-nearest neighbor coupling (the J1J2J3 model), making use of 
Takahashi's modified spin-wave (MSW) theory supplemented by ordering vector op- 
timization. We compare the MSW results with exact diagonalization and projected- 
entangled-pair-states calculations, demonstrating their qualitative and quantitative 
reliability. We find that MSW theory correctly accounts for strong quantum effects on 
the ordering vector of the magnetic phases of the models under investigation: in par- 
ticular collinear magnetic order is promoted at the expenses of non-collinear (spiral) 
order, and several spiral states which are stable at the classical level, disappear from 
the quantum phase diagram. Moreover, collinear states and non-collinear ones are 
never connected continuously, but they are separated by parameter regions in which 
MSW breaks down, signaling the possible appearance of a non-magnetic ground state. 
In the case of the spatially anisotropic triangular lattice, a large breakdown region 
appears also for weak couplings between the chains composing the lattice, suggesting 
the possible occurrence of a large non-magnetic region continuously connected with 
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the spin-liquid state of the uncoupled chains. 



PACS numbers: 75.30.Ds,75.30.Kz,75.10.Jm,75.50.Ee 



I. 



INTRODUCTION 



Low- dimensional frustrated quantum spin systems can display an intriguing interplay 
between order and disorder: classical order has been shown to be quite resilient in two or 
three dimensions [IHlj; frustration, however, can lead to the melting of magnetic long-range 
order (LRO) and the emergence of quantum disordered states like valence-bond solids or 
resonating valence bond states [5j [6]. Understanding such magnetically disordered quan- 
tum phases is important for the search for fractionalized excitations in two dimensions [5j, 
as well as for the understanding of the behavior of layered magnetic insulators/metals in 
which magnetism is disrupted by charge doping, leading to dramatic phenomena such as 
superconductivity at high critical temperature [7H9]. 

A large variety of magnetic materials can be described by the Heisenberg Hamiltonian 



where Si is a quantum spin-5 operator at site i. In this paper, we will focus on the 
antiferromagnetic case for 5 = 1/2, and on two-dimensional frustrated lattices. Quasi- 
two-dimensional frustrated antiferromagnetism is relevant to a variety of 5 = 1/2 com- 
pounds, realizing the spatially anisotropic triangular lattice (e.g., in CS2CUCI4 [TU] and 
/t-(BEDT-TTF)2Cu2(CN)3 [TT| H~2] . etc.), or the frustrated (J1J2) square lattice (e.g., in 
Li 2 VOSi(Ge)0 4 , VOM0O4 [13J, BaCdVO(P0 4 ) 2 [14J, etc.). For both lattice geometries, the 
Heisenberg model is expected to display spin-liquid phases for particular values of the frus- 
trated couplings, although the extent and nature of these spin-liquid phases is still under 
theoretical debate, both for the spatially anisotropic triangular lattice (SATL) [T0H2T] and 
for the frustrated square lattice j22H2H]- 

In this work, we investigate the 5 = 1/2 Heisenberg antiferromagnetic Hamiltonian on 
two-dimensional frustrated lattices making use of Takahashi's modified spin-wave (MSW) 
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theory |29j . supplemented with the optimization of the ordering vector |30j. In a previous 
paper |31j . we have shown that (for the SATL with XY interactions) this approach provides 
a significant improvement over conventional spin-wave theory (as well as over conventional 
MSW theory), as it allows to correctly account for the dramatic quantum effects occurring 
to the form of order which appears in frustrated quantum antiferromagnets, and for the 
quantum corrections to the stiffness of the ordered phase. In particular, a very low stiffness, 
or the complete breakdown of the theory, provide strong signals that the true ground state 
might be quantum disordered; hence, this method serves as a viable approach to finding 
candidate models potentially displaying spin-liquid behavior. For a more detailed description 
of the formalism we refer the reader to Ref . |3T] . 

Here, we apply this MSW theory with ordering vector optimization to the Heisenberg 
SATL, as well as to the square lattice with nearest, next-to-nearest and next-to- next-to- 
nearest neighbor couplings (the J1J2J3 model |27| I32H31] ) . Both models feature a very 
complex T = phase diagram, with spirally and collinearly ordered regions, whose ordering 
vector is subject to strong quantum corrections with respect to the classical (S — > 00) limit. 
They also feature extended breakdown regions for MSW theory, pointing at the possible spin- 
liquid nature of the true ground state of the system. Comparison with numerical results 
coming from exact diagonalization and projected-entangled-pair-state (PEPS) calculations 
show that MSW theory correctly accounts for some of the most salient features of the 
quantum phase diagram of these systems, and that it hence represents a very versatile tool 
to probe the robustness (or the breakdown) of a semi-classical description of the ground 
state of frustrated quantum magnets. 

The remainder of this paper is organized as follows: Section [IT] presents the ground state 



phase diagram of the SATL with nearest-neighbor Heisenberg interactions; in Section III 



we calculate the ground state phase diagram of the J1J2J3 model; finally, in Section IV 



we present our conclusions. The technical aspects of MSW theory applied to Heisenberg 
antiferromagnets are presented in the Appendix. 
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II. MSW THEORY ON THE SPATIALLY ANISOTROPIC TRIANGULAR 
LATTICE WITH NEAREST-NEIGHBOR HEISENBERG-BONDS 

The triangular lattice with Heisenberg interactions has been considered as one of the 
first candidate systems for quantum-disordered behavior in the ground state [5]. Recently, 
the phase diagram of the spatially anisotropic triangular lattice (SATL) up to values of 
a = t 2 /ti = 1 has been studied by Yunoki and Sorella using variational quantum Monte 
Carlo methods [16J. They find that the gapless spin-liquid phase of the isolated chains 
(t 2 = 0) persists also at finite coupling up to a critical value a ~ 0.65, followed by a gapped 
spin liquid; for a ~ 0.8 the gap closes and the system undergoes an ordering transition to 
spiral order, continuously connected with the 3-sublattice order of the isotropic Heisenberg 
antiferromagnet (a = 1). This scenario is still controversial, however: studies based on 
low-energy effective field theory for the description of the coupled chains in the case a < 1 
indicate that the system might still exhibit long-range antiferromagnetic order even for very 
weak coupling among the chains. This form of order results from high-order perturbation 
theory in the inter-chain coupling, and it is necessarily very weak, given that numerical 
methods cannot detect it. Its observation is clearly beyond the capabilities of our MSW 
approach. Coming from the large-a limit, series expansions by Weihong et al. indicate 
that 2D-Neel order - appearing on the square lattice defined by the dominant /^-couplings 
- persists down to a ~ 1.43, followed by a phase without magnetic order in the interval 
1-1 ~ a ~ 1-43 |15| . Below this region the authors find incommensurate spiral order 
connecting continuously to the isotropic point a = 1. In Ref. [35j, qualitative similar results 
have been obtained using the Schwinger-boson approach. The resulting phase diagram 
differs strongly from the classical one, which is characterized by spiral order for < a < 2, 
and by Neel order for a > 2. The classical phase diagram is contrasted with the quantum 
mechanical one (composed from Refs. jT3] and |16j ) in Fig. [T] It is interesting to notice that 
a qualitatively similar phase diagram has been obtained recently by some of us for the XY 
model on the SATL pTJ, [36]. 

A variety of experiments have been carried out on magnetic compounds described by the 
Heisenberg model on the SATL, with results that are still controversial. For instance neutron 
scattering experiments of Coldea and coworkers |TU] on CS2CUCI4, where a ~ 1/3, claimed 
evidence that the low-energy physics is governed by spinons, fractionalized excitations with 
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~ 0.6 ~ 0.8J ~ 1.1 ~ 1.43 tijtx 

spiral order 



Figure 1: (a) Classical ground state phase diagram of the SATL with a sketch of the ID state at 
a = 0, the spiral state at a = 1, and the 2D-Neel state for a > 2. (The horizontal black bonds 
have strength t±, and the diagonal red bonds have strength £2-) (b) The quantum mechanical phase 
diagram changes considerably due to order-by-disorder effects and the appearance of spin liquids 

PUCE]. 

S = 1/2 which represent the elementary excitations in the case of uncoupled chains. Yet, 
Ref. [19] showed that, for a finite inter-chain coupling, spinons tunnel between chains in 
bound pairs with S — 1 (so-called triplons), so that the fractionalization in two dimensions is 
strictly speaking not present. Ref. [19] argues that the spinons in CS2CUCI4 are descendants 
of the excitations of the individual ID chains and not characteristic of any exotic 2D state. 
This further reinforces the idea of a quasi one-dimensional behavior up to relatively high 
inter-chain interactions mentioned in the previous paragraph. 

A. MSW predictions for the ground-state phase diagram 

In this section, we discuss the ground-state phase diagram resulting from the predictions 
of MSW theory for the 5=1/2 SATL with nearest-neighbor (NN) Heisenberg interactions. 

In order to assess the validity of MSW results, we compare them with exact diagonaliza- 
tions (ED). Using the Lanczos method, we compute the ground state of small clusters of 14, 
24, and 30 spins. The considered geometry for the 30-spin system can be found in Fig. [2} 
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Figure 2: Cluster of 30 spins for which we carried out ED. The 24-spin system is equivalent, only 
with the top and bottom rows removed. Black dots denote sites, the horizontal black bonds have 
strength t\, and the diagonal red bonds have strength ti- 

The 24-spin system can be obtained from it by removing the top and bottom rows. The 
14-spin cluster is an equivalent system with rows of 2, 3, 4, 3, and 2 spins. The clusters are 
chosen for their symmetry with respect to reflection along the coordinate axis, and for their 
ratio of the number of ^-bonds (red) to the number of tx-bonds (black), which lies close 
to the bulk value of 2. We use open boundary conditions to allow for the accomodation of 
spiral order with arbitrary wave vector. 

We find that, due to the peculiar geometries chosen, there exist parameter ranges where 
the ground state falls into the threefold degenerate triplet with total spin equal to unity. 
Nonetheless, we restrict our calculations to the M* otal = subspace (with M* otal being the 
z component of the total spin), and the M* otal = ±1 states are excluded. This results in an 
apparent breaking of the x—z symmetry (the x-y symmetry is preserved). This symmetry 
would be recovered by averaging over the whole triplet subspace. The reason for such an 
apparent symmetry breaking resides in the particular geometry of the cluster considered, 
which complicates the comparison between different system sizes. This triplet physics might 
play an important role for bigger systems, although one cannot draw conclusions about the 
thermodynamic limit from the small clusters considered. A non-trivial triplet physics could 
be especially an issue for variational methods restricting their focus to the singlet subspace. 

The lattice sizes considered in the MSW calculations are 32 x 32 spins and the infinite 
lattice limit, which is achieved by transforming sums over the first Brillouin zone into inte- 
grals. Figures [3] to [8] show that the difference between the lattice sizes is insignificant except 
near quantum phase transitions, which is expected because of the divergence of correlation 
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lengths near criticality. 

1. Parameter regions where MSW theory fails to converge. 
Convergence in the self-consistent equations of MSW theory with ordering vector opti- 



mization, Eqs. (A3- A7 A9), cannot be achieved in selected regions of the ground state phase 
diagram, namely for a < 0.65 and for 1.14 < a < 1.3. (Interestingly, convergence is restored 
in the pure ID limit, a = 0, for which the theory formulates surprisingly good predictions.) 
This breakdown of convergence corresponds to the appearance of an imaginary part in the 



spin-wave frequencies, Eq. (A5), signaling an instability of the ordered ground state. The 
breakdown of a self-consistent description of the system in terms of an ordered ground state 
is strongly suggestive of the presence of a quantum-disordered ground state in the exact 
behavior of the system. Hence, one can interpret these parameter regions as candidates for 
the spin-liquids predicted from Refs. [T5l[T6] [compare Fig. |l](b)]. Both for a < 1 and a > 1, 
we find that the breakdown region of MSW appears to be fully contained within the region 
of SL behavior (either gapped or gapless) estimated in Refs. [TBI ITS] . Hence MSW theory 
is seen to possibly underestimate the width of the quantum-disordered regions in the phase 
diagram, which is to be expected due to the partial account of quantum fluctuations given 
by MSW theory. 

2. MSW ground state energy in comparison with previous results. 

Table |TJ demonstrates that the energy from MSW theory compares very well to results 
that were obtained recently by Yunoki and Sorella by variational quantum Monte Carlo 
methods [16] . also plotted in Fig. |3} For comparison, we also show the curve that they 
obtain with a projected-BCS (p-BCS) wave-function. In the isotropic triangular lattice, 
the MSW energy compares also favorably to the data from the Green's function Monte 
Carlo method with stochastic reconfiguration (GFMCSR) from Ref. [10], but both energy 



and order parameter (see section II A 3 ) lie closest to the variational quantum Monte Carlo 
calculation from Weber et al. [38J, who used a mixture of a BCS wave-function and a wave 
function with spiral order as their starting point (BCS+spiral). 

At a = the MSW value Eq (a = 0) = —0.4647 is relatively close to the exact result 



8 



Method 


a = 


a = 0.7 


a = 0.8 


o = l 


a = 00 


M at a = 1 


exact, thermodynamic limit 


-0.443147 












exact, N = 30 (present study) 


-0.4127 






-0.5471 




0.1314 


exact, N = 40, extrapolated |37j 










-0.6701 




VMC (RVB) [16] 








-0.5123(1) 




0.0 


VMC (RVB with u = 0) ITB1 








-0.5291(1) 




0.0 


VMC (BCS+spiral) [38] 








-0.532(1) 




0.36 


VMC (p-BCS) Hi 


-0.442991 


-0.46467 


-0.47840 


-0.5357(1) 




0.0 


FN [16] 




-0.47051 


-0.48521 


-0.53989(3) 




0.1625(30) 


FNE [16] 




-0.47171 


-0.48691 


-0.54187(6) 




0.1765(35) 


GFMCSR [39] HQ] 








-0.545(2) 




0.205(10) 


series expansion 1391 a 










-0.6696(3) 




LSW H6][39] a 








-0.538(2) 




0.2387 


MSW (present study)" 


-0.4647 


-0.4639 


-0.4775 


-0.5303 


-0.6699 


0.3426 



These methods 



do not provide a rigorous upper bound for the ground state energy. 



Table I: Comparison of the ground state energy per spin derived by various methods, for some 
values of a. VMC stands for 'variational quantum Monte Carlo' where the wave function used is 
given in brackets |16| [38] . FN is short for lattice fixed node and FNE for the improved FN effective 
Hamiltonian method [16j. Furthermore included are the exact result of the Heisenberg chain in the 
thermodynamic limit, the ED results for the 30-spin cluster, and the ED results from finite size 
extrapolations of calculations on clusters of up to 40 spins [37J. Also given are the estimates of LSW 
theory from Ref. [16] and the Green's function Monte Carlo method with stochastic reconfiguration 
(GFMCSR) [40] , The last column gives the staggered magnetization or, in the case of MSW theory, 
the population of the zero mode Mq. 



of the one-dimensional case, —(In 2 — 1/4) = —0.44315. However, it is located below the 
exact value. This apparent puzzle is easily resolved by noticing that the MSW method is 
not variational due to the incomplete inclusion of the kinematic constraint (see Appendix). 
We also notice that the ground state energies derived from ED of the 30-site system lie very 
close to the values from the other methods except in the ID phase. This could be attributed 
to the small system size: if the interpretation is correct that for small a the Heisenberg 
SATL is in a ID-like phase with algebraic correlations, it is natural that finite size effects 
play a very important role in the critical ID phase. This would explain the strong deviation 
of the ED energy in that parameter region. 

On the square lattice (a — > oo), Takahashi showed already twenty years ago the extremely 
good performance of MSW theory j29]: its ground state energy per spin is —0.6699, which 
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Figure 3: MSW results for the ground state energy lie close to results from previous studies. 
Shown are the data of Ref. |16j for their variational quantum Monte Carlo (VMC) Ansatz with a 
projected BCS wave-function (p-BCS) and the improved FN effective Hamiltonian method (FNE). 
We further display the value obtained in the isotropic limit by Ref. [38] by use of a VMC method 
with a mixture of a BCS and a spiral ordered wave- function (BCS+spiral), and the exact result of 
the ID limit. The numbers in the labels of the curves are the respective system sizes considered. 

is in excellent agreement with the QMC result —0.669437(5) |41| . 



3. Order parameter and spin stiffness from MSW theory. 

Our next step is to determine the regions where the presence of a finite order parameter 
Mo and spin stiffness T reveal magnetic long-range order (LRO). Even when Mq and T are 
finite, a caveat is still in order: a finite order parameter with a very small stiffness might 
suggest that taking quantum fluctuations more completely into account than in MSW theory 
could lead to a completely disordered state. 

The order parameter Mq, drawn in Fig. |4j shows that magnetic LRO is present in the 
intervals 0.65 < a < 1.14 and a > 1.3. This is to be contrasted with linear SW (LSW) 
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Figure 4: M$ from MSW theory compared with the LSW values and ED results (see section II A 5). 



The numbers in the labels of the curves are the respective system sizes considered in the calculations. 

theory, which predicts the breakdown of magnetic order only for a < 0.3 [42J. However, in 
the isotropic case, a — 1, MSW theory predicts a stronger order parameter than what is 
predicted by LSW, as well as by most of the other numerical estimates, which are presented in 
Table [Tj In the square lattice limit, a — > oo, on the other hand, both MSW and LSW theory 
attain a staggered magnetization of 0.303, which compares favorably with the most recent 
estimates M = 0.311, based upon diagonalizations of small clusters of up to 40 spins [37J. 
The MSW order parameter drops drastically upon approaching the region 1.14 < a < 1.3 
and when reaching the region a < 0.65 from above, the regions where the self-consistent 
description breaks down, further corroborating the assumption that in these regions magnetic 
LRO disappears in the true quantum ground state. This assumption is strongly reinforced 
by considering the Gaussian spin stiffness (Fig. [5]): It vanishes at a = 0.65 and it drops 
significantly when approaching a = 1.14 from below. 

There are various special cases of the SATL for which the spin stiffness has been calculated 
previously. In the square lattice limit, a — > oo, MSW theory gives p\\/a = 0.216, somewhat 
overestimating the value from QMC p\\/a = 0.175(2) [41J. In the isotropic triangular lattice, 
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Figure 5: (a) Gaussian spin stiffness T (for the 32 x 32 and the infinite system) and (b) the 
components of the spin stiffness tensor (for the infinite system). The mixed second derivative 
p xy vanishes for symmetry reasons. The curves labeled 'partial' were obtained by application of 



Eq. (A10). 



a = 1, the spin stiffness from the MSW approach is p\\/a = 0.113. This value falls between 
the LSW spin stiffness p\\/a = 0.122 (Ref. jl3]) and the estimate obtained from ED calcu- 
lations after finite size extrapolation, p\\/ot = 0.075 [43j. In the limit of decoupled chains, 
a = 0, MSW theory achieves convergence (which was lost in the interval < a < 0.65) and 
provides a spin stiffness p xx /a = 0.309 in the thermodynamic limit, relatively close to the 
exact result in the thermodynamic limit, p xx /a = 1/4 p 
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4- Spin and chirality correlations from MSW theory 

Now we describe the ordered phases found by the MSW Ansatz for the Heisenberg SATL 
in more detail. To this end, we analyse the following quantities 



1. The ordering vector Q (Fig. [6J). Three limiting values for the ordering vector are 
known. For a = intra-chain antiferromagnetic (Neel) order is described by Q = net. 
For a — > oo square-lattice Neel order is described by Q = 2ttx. In the isotropic lattice 
(a = 1), the threefold symmetry forces the ordering vector to Q = ^-x; 



2. The spin-spin correlations (Fig. [7]). We analyze the spin-spin correlations of nearest 
neighbors through the two-site total spin, 

T ij = ±({S i + S j ) 3 ) = (S i -S j ) + l. (2) 

This quantity vanishes if the spins are in a singlet, which is equivalent to perfect 
anticorrelation, takes the value | if they are uncorrelated, and the value 1 if the spins 
form a triplet, which means perfect correlation; 



3. The mean chiral correlations (Fig. [8]). Spiral phases carry not only a magnetic 
order parameter, but also a chiral order parameter. In particular, a vector chi- 
rality can be defined on an upwards pointing triangle with counter-clockwise la- 
beled corners k) as 03] ka = -^m [Si X Sj + Sj x Sk + Sk x Si] z , and on a 
downwards pointing triangle with counter-clockwise labeled corners as Ky = 

[Si x Si + Si x Sj + Sj x Si] z . Chirality correlations are defined as [IB] 

4>- = ({k-a - «v) {ka< - «v)) > ( 3 ) 



where the triangle pairs (A, V) and (A', V') share a T\ = (1,0) edge. In Fig. [8j we 
plot the average chirality correlation of the central plaquette with all other plaquettes, 
normalized to the theoretical maximum 4/9. The MSW data have been obtained by 
expanding the chiral correlation up to the fourth order in the boson operators, which 



is consistent with the truncation of the bosonic Hamiltonian Eq. (Al) to the same 
order. Going to higher orders does not change the outcome in the regions where M 
is large, but can yield different results where M is small. 



13 



2 

%1.5 



0.5 1 1.5 2 

Figure 6: First component of the ordering wave- vector, Q x , from MSW theory. Also shown are 
the classical values and the ED results. The black circle marks the order vector Q x = 120° of the 
isotropic triangular lattice which is attained classically and by the spin- wave theories at a = 1. 
The numbers in the labels give the system sizes. 

A comparison of these quantities shows a spiral phase at around 0.65 < a < 1.14 and 
a 2D-Neel ordered phase for a > 1.3. Moreover, when approaching a ~ 0.65 from above, 
the ordering vector, the spin-spin correlations and the ground state energy approach their 
respective ID values. This is an indication that below a ~ 0.65 the true ground state of 
the system may enter a ID-like spin-liquid phase. Nonetheless, the vanishing of the spin 
stiffness for a — > 0.65 + is not consistent with the onset of a gapless ID spin-liquid phase, 
for which the spin stiffness should remain finite. Hence, the MSW results rather suggest 
that the phase appearing below a = 0.65 is a gapped spin liquid, and that the gapless ID 
spin-liquid phase, connected continuously with the limit a = 0, is only attained for even 
smaller a. This seems consistent with the prediction of Ref. [16] that a gapped spin-liquid 
phase separates the spirally ordered phase from the ID-like gapless disordered one. 

5. Order parameter and correlations in comparison with exact diagonalization. 
In the case of ED, the static structure factor 

sa = E ( s ? s ?) e ' tk ' r ' J (« = y> z ) ( 4 ) 
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Figure 7: MSW and ED results show similar behavior of Jo jTi , with Tj = T\ = (1,0) (solid lines) 
and Tj = T2 = (l/2, \/3/2) (dashed lines), respectively. The numbers in the labels give the system 



sizes. 
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Figure 8: Comparison of the MSW and ED results for the mean chiral correlation normalized to 
the theoretical maximum of 4/9. The numbers in the labels give the system sizes. 



allows to extract the order parameter M a , which is defined as M a = \JS a (Q), where Q 
is the ordering vector associated with a peak in S a (k). In the thermodynamic limit, this 
is the equivalent to Mo from MSW theory. A comparison of both quantities can be found 
in Fig. [4} We plot both M x and M z due to the anisotropy caused by the triplet physics 
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mentioned at the beginning of this section. Discontinuous jumps in the ED magnetizations 
are due to the change of the spin sector hosting the ground state, going from the singlet 
sector (characterized by M x = M z ) to the triplet sector (characterized by M x ^ M z ). We 
observe very severe deviations between the ED data on the one side and the predictions from 
LSW and MSW theory on the other side: in particular, apart from the deviations between 
M x and M z , the ED data appears to be almost constant over a large a interval. The strong 
difference between ED results on the one hand and MSW/LSW predictions on the other 
can also be attributed to very significant finite-size corrections to the ED data - finite-size 
effects are particularly pronounced here, due to the open boundary conditions of ED clusters. 
Nonetheless, for a = 1 the magnetization of the 30-site cluster gives M x = M z « 0.13, lying 
close to recent Monte Carlo estimates [21J. 

From the location of the peak of the structure factor one can extract the vector of pre- 
dominant ordering, Q, the x-component of which is plotted in Fig. [6] Remarkably, for the 
30-site cluster the Q corresponding to M x (labeled as Q x in the figure) indicates a transition 
from spiral to Neel order at around a w 1.4, which lies well below the classical threshold 
a = 2. On the contrary, the Q corresponding to M z (labeled as Q z ) increases smoothly up 
to a ~ 2, where it undergoes a discontinuous transition to the square-lattice Neel value as 
well. However, increasing the system size from 24 to 30 spins shifts significantly the curves 
of Q x and Q z to the left, suggesting that for even larger sizes both curves might exhibit a 
discontinuous transition to the Neel ordering vector for a value of a close to the transition 
indicated by MSW, a ~ 1.3. Finally, we notice that at a = 1 the ED results deviate from 
the isotropic value Q x = 120° because the required threefold symmetry is broken by the 
shape of the simulation cluster, Fig. [2j 

The nearest-neighbor spin-spin correlations T y -, Eq. ([2]), [71 J are in qualitative agreement 
with the MSW results as well (Fig. [7]). In particular, they show ID-like behavior at small 
a, a spiral phase in an intermediate parameter range around the isotropic limit a — 1, and 
a 2D-Neel structure at large a. 

Finally, we focus on the chirality correlations. Comparing such correlations for the 14, 24, 
and 30 spin clusters shows that they are strongly suppressed for a < 0.5 and for a > 1.4 when 
going to larger lattice sites. This indicates that a non-spiral phase appears in this region 
in the thermodynamic limit, in agreement with our MSW calculations. The persistance of 
significant correlations in the region 0.5 < a < 1.4 indicates that spiral order in the ground 
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state might persist in a portion of this parameter range. 

In summary, despite the significant deviations in the magnitude of the order parameter, 
both ED and MSW theory give a coherent picture, both qualitatively and quantitatively, of 
the evolution of the nature of spin-spin correlations upon increasing the a parameter, going 
from quasi-lD to spiral to Neel. 

B. Discussion 

Despite its limitations, the MSW approach with ordering vector optimization reproduces 
faithfully the main characteristics of the phase diagram as sketched in Fig. [T] (b), and 
thus remarkably improves on the results that were previously obtained for this model with 
conventional spin-wave theories. A breakdown of magnetic order - along with a variety of 
observables like the ordering vector or nearest-neighbor spin-spin correlations - indicates 
that a ID-like spin liquid might be attained below a ~ 0.65. Due to the partial account of 
quantum fluctuations provided by MSW theory, we can safely take this as a lower bound for 
a spin liquid in the true ground state. Furthermore, we find a relatively small region with 
spiral LRO between 0.65 < a < 1.14. For a > 1.30 the system is ordered at the 2D-Neel 
wave- vector. Between 1.14 < a < 1.30 the breakdown of convergence suggests another 
candidate region for spin-liquid behavior. 

III. MSW THEORY ON THE Jihh MODEL 

In this section, we investigate another paradigmatic frustrated spin model, the J1J2J3 
model on the square lattice. It involves couplings between nearest- neighbors (NN), Ji, next- 
nearest- neighbors (NNN), J2, and next-next-nearest-neighbors (NNNN), J3. A sketch of the 
geometry of the system may be found in Fig. [9] (a). This model allows to continuously tune 
the Hamiltonian from an unfrustrated antiferromagnetic square lattice to a highly frustrated 
magnet. 
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Figure 9: (a) A detail of the geometry of the J1J2J3 model on a square lattice. Nearest neighbors 
are coupled with bonds of strength J\ (black), next-nearest neighbors (along the diagonals) with J2 
(blue) and next-next- nearest neighbors with J3 (red), (b) The classical phase diagram of the J1J2J3 
model shows four ordered phases: Phase I is characterized by Neel order on the square lattice. In 
phase II the system decouples into two independently Neel ordered sublattices with a doubled unit 
cell each. Phases III and IV are spirally ordered with Q = (q,ir) and Q = (q,q), respectively. 

A. Classical and quantum mechanical phase diagram of the J1J2J3 model at T = 



The classical phase diagram of the J1J2J3 model [Ml B71 - H9] is sketched in Fig. [9] (b). 
One identifies: 

I) A 2D-Neel phase with Q = (jr, tt) just as in the unfrustrated square lattice. It is 
delimited by the classical critical line (J 2 + 2J 3 ) / J\ = 1/2; 

II) A phase where the system decouples into two independent J2— sublattices with a dou- 
bled unit cell. Both are Neel ordered individually. This phase is infinitely degenerate 
because the two sublattices can be rotated one with respect to the other without 
affecting the energy; 

III) A spiral phase with ordering vector Q = (q, 7r), where q varies continuously over the 
phase diagram; 

IV) A second spiral phase, this time with ordering vector Q = (q,q); q — >■ it/2 for J3 — > 00, 
attaining the limit of two decoupled and Neel-ordered J 3 — sublattices. 

This phase diagram is believed to change considerably in the quantum limit [271 132HM] : 
In phase II quantum fluctuations select the columnar ordered states with Q = (tt, 0) or 



18 



Q = (0, tt) from all the possible classical states. Furthermore, the Neel phase I increases in 
size considerably and Neel order persists up to the vicinity of the line (J2 + ^3) I J\ — 1/2. 
In the vicinity of this line, the classical order is believed to be destabilized and to be replaced 
by a non-magnetic state. The controversy about the exact nature of the ground state in this 
highly frustrated region, however, is still not settled. In particular, it has been suggested 
that it could have the nature of a columnar valence bond crystal [50j with both translational 
and rotational broken symmetries, of a plaquette state with no broken rotational symmetry 
|27j, or of a spin liquid with all symmetries restored [o"T1 - I55| . 

In the following, we investigate the quantum model using the modified spin-wave (MSW) 
formalism, and compare it to recent results from projected entangled-pairs states (PEPS) 
calculations. The MSW lattice size is again N = 32 x 32. In most of parameter space, a 
lattice of N = 32 x 32 spins is essentially already converged to the infinite lattice, except 
close to a quantum critical point. 

In Ref. |56j . some of us reported numerical calculations of the J1J2J3 model based on the 
projected entangled-pair state (PEPS) variational Ansatz for varying lattice sizes. In the 
following, we will focus on the extrapolations to the thermodynamic limit, except if stated 
otherwise. 

We first discuss in more detail the special cases of the J1J2 model (i.e., J3 = 0) and the 
J1J3 model (i.e., J2 = 0). Both models have been studied before within the MSW formalism 
[3T)| loTHoT] . On the one hand, we confirm existing results on the J\ J 2 case, for which the 
optimization of the ordering wave-vector returns only two possible values (corresponding to 
Neel order [Q = (tt, tt)] or columnar order [Q = (tt, 0) or Q = (0,7r)]), and we give further 
insight into the spin stiffness and the dimer-dimer correlation functions. On the other hand, 
we analyze the J\ J3 model with optimization of the ordering wavevector, which proves crucial 
to correctly capture the quantum effects on the classical spiraling phases appearing in this 
case [30J . Finally, we give an overview of the entire quantum ground state phase diagram of 
the J1J2J3 model. 



B. Ground state properties of the J\ J2 model 



Figures 10 and 11 report the results for the J1J2 model from the MSW method as well 
as from PEPS calculations. For comparison, we also plot the values for the energy and 
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Figure 10: For the J\ J2 model the y-component of the ordering vector shows a considerable shift 
in the quantum model with respect to the classical value. 

magnetization that where obtained in Ref. [37] from diagonalization of small clusters. In 
agreement with other methods, e.g., exact diagonalization (ED) [371 1621 [63] or Schwinger 
bosons [64J, MSW theory finds Neel order with Q = (7T,7r) at small J2/J1 and columnar 



order with Q = (7r,0) or Q = (0, 7r) at large J2/J1 (see Fig. 10). As it is well known 
from previous studies, there is a region between 0.56 < J2/J1 ^ 0.62 where the 2D-Neel 
ordered and the columnar state are both stable solutions within MSW theory. The starting 
point of the self-consistent calculations determines which type of order is returned as the 
solution. However, the solutions differ in energy and therefore one of them is only a local 
free energy minimum of the self-consistent equations. The transition from 2D-Neel order to 
columnar order takes place at J2/J1 — 0.6. For the PEPS results, we extract the wave vector 
of dominant spin correlations Q PEPS from the location of the peak of the static structure 
factor, 

M(k) = 




^J^iSi-Sje-^). (5) 



In agreement with the MSW prediction, Q PEPS is located at the Neel value (7r,7r) up to 
J2/J1 = 0.6, while above this it lies at the value of columnar order (7r, 0). 

We find a remarkable correspondence of the ground state energy per spin between 
the MSW prediction and ED results extrapolated to the infinite lattice from Ref. [37] 



[Fig. 11 (a)]. Moreover, the noticeable kink associated with the Neel-to-columnar transi- 
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Figure 11: Comparison, for the J% J 2 model, of MSW data to PEPS results extrapolated to the 
thermodynamic limit and ED from Ref. |37j [for the 40-spin cluster (labeled 'exact 40') and extrap- 
olated to the thermodynamic limit ('exact 00')]. For the MSW method, the curves obtained when 
starting the self-consistent iteration from a Neel state (thick red line) and from a columnar ordered 
state (thick dot-dashed green line) are both included. The figures show (a) the ground state energy 
of the central spin; (b) the MSW order parameter Mq compared to M (tt, tt) (Neel) and M (tt, 0) 
(columnar) derived from PEPS calculations; (c) the Gaussian spin stiffness, and (d) components 
of the spin stiffness tensor. In the Neel phase p xx = p yy by symmetry. The partial spin stiffnesses 
p^ rtial are found to equal the total ones, p a @. 



tion of MSW theory at J 2 / J\ = 0.6 is exhibited as well by the 40-sites system from Ref. 
Therefore, ED confirms that J2/J1 = 0.6 marks a transition point, although in the true 
ground state such a transition might connect the columnar state to a quantum-disordered 
state. A similarly good agreement is found with the PEPS results extrapolated to the infinite 
size limit. 



As shown in Fig. 11 (b), at small J2/J1, i-e-, deep in the Neel phase, the finite size extrap- 
olation of the ED staggered magnetization from Ref. [37J lies very close to the MSW results. 



As it is well known |29j . in the unfrustrated square lattice limit (J 2 = 0) the MSW value 
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Mq = 0.303 is only slightly smaller than M = 0.311 from ED. For the PEPS calculations an 



analogous quantity can - similar to section |II A5| - be derived from the peak height of the 



static structure factor, Eq. (pi). We show its finite size extrapolation in Fig. 11 (b). In the 
Neel phase PEPS agrees very well with MSW theory, considerably better than ED, which 
decreases faster towards the strongly frustrated region. In the entire columnar phase, how- 
ever, PEPS and ED data lie closer together, while MSW overestimates the order parameter. 
Around the transition, however, agreement between PEPS and MSW theory is very good. 
The PEPS data suggest that the magnetically disordered region, predicted by ED to occur 
in the range 0.35 < J2/J1 ^ 0.66, is either much smaller or does not occur at all. 

The MSW spin stiffness p\\ = (p xx + p yy ) /2, however, while being finite for any considered 
value of the ratio J2/J1, is strongly suppressed in the region 0.3 < J2/J1 ^ 0.6 [Figs. 11 (c) 
and (d)], suggesting as usual that accounting for quantum fluctuations beyond the MSW 
approximation could lead to the disappearance of magnetic order. A suppression of spin 
stiffness is also observed in previous results coming from ED of finite clusters [62J or from 
the Schwinger boson approach |6T| [55] . As a consequence, even though MSW admits a 
stable solution with magnetic order for any J2/J1 value, for J2/J1 = 0.6 it exhibits a clear 
transition from soft Neel order to a stiff columnar order, suggesting that this transition 
could actually separate the columnar state from a quantum disordered phase. 

1. Dimer correlations in the J1J2 model. 

The nature of the state in the transition region between Neel and columnar order, where 
magnetic order is strongly reduced, can be further investigated through the study of the 
dimer-dimer correlations 

C ijkl = ((# • Sj) (S k ■ S t )) , (6) 

where k and I, and i and j are pairs of neighboring spins. Figure |T2| sketches the expectation 
for the dimer-dimer correlations in (a) a columnar valence bond crystal and (b) a columnar 
magnetic state. 



In Fig. 13, we show the spatially resolved dimer-dimer correlations from MSW theory. 
Below J 2 / Ji = 0.6 the dimer-dimer correlations have a structure compatible with a Neel state 
(namely they are positive and nearly equal for all bond pairs), while above J2/J1 = 0.6 the 
dimer-dimer correlations acquire the expected structure in a columnar state, with opposite 
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Figure 12: Sketch of dimer-dimer correlations in (a) a valence bond crystal and (b) a columnar state. 
Black dots are lattice sites. Blue solid (red dashed) lines are dimers correlated (anti-correlated) 
with the central dimer (thick black line). 

signs for the correlations between dimers of the same spatial orientation (both horizontal 
and both vertical) and between dimers of opposite orientations. Nonetheless, for J%/ J\ < 
0.7, remarkably MSW theory shows a short-range modulation in the strength of the dimer 
correlations whose structure is compatible with that of a valence bond crystal. Although 
MSW theory is not appropriate to characterize non-magnetic states such as a valence bond 
crystal, it is remarkable to observe that it identifies a columnar valence-bond structure as 
the dominant form of dimer correlations at short range; this indication is consistent with, 
e.g., the results of PEPS [66j, which also point towards columnar valence-bond order in the 
non-magnetic region of the J1J2 model. 



C. Ground state properties of the J1J3 model 

We now turn to the J\Jz model. Classically, this model has a transition from Neel to 
spiral order at J 3 = 0.25Ji. Recent PEPS calculations show that for S = 1/2 Neel order 
persists up to approximately J3/J1 = 0.3 [56]. Above this point the peak of the structure 
factor is still at the Neel ordering vector (tc, it) but its height vanishes in the thermodynamic 
limit, which suggests a complete loss of magnetic LRO. A different type of LRO arises anew 
at approximately J3/J1 = 0.6 with an ordering vector Q = (q,q) that tends to (7r/2,7r/2) 



in the limit of large J 3 (see Fig. 14). For large enough J 3 the nature of the ordered phase 
becomes similar to that of the classical limit. 

The optimization of the ordering wave- vector within MSW calculations shows that, for 



small J3/J1, Neel order persists up to J3/J1 = 0.39 (see Fig. 14), confirming the assump- 
tion that quantum fluctuations stabilize Neel order against spiral order with respect to the 
classical limit. Coming from the opposite limit of J3 ~ Ji, we observe a spiral phase with 
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Figure 13: MSW correlations of the black central dimer with the other dimers of a 32x32 lattice 
(zoom on central region) . The thickness of the lines is a non- linear function of the absolute strength 
of the dimer correlations. Note the change of the maximum of the linear color-scales for different 
values of J%l J\. Below Jij J\ = 0.4 and above Jij J\ = 0.9, the qualitative changes are minimal. 

continuously varying pitch vector Q = (q,q), where q approaches 7r/2 for J3/J1 — > 00, and 
increases up to q w 0.77T for J3/J1 — >■ 0.52 + . In the region 0.39 < J3/J1 < 0.52, conver- 
gence of the MSW calculations breaks down, which points at a possible spin-liquid phase, 
in agreement with the predictions from PEPS calculations. 



Fig. 15 (a) shows the PEPS energy extrapolated to the thermodynamic limit. Agreement 



to the MSW results is again found to be extremely good. 
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Figure 14: Position (Q x , Qy) of the peak of the structure factor for PEPS and the ordering vector 
Qx = Qy of MSW theory for the J1J3 model. A comparison to the classical ordering vector 
Q x l = shows that quantum fluctuations stabilize Neel order. 
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Figure 15: MSW and PEPS results on the J1J3 model. Shown is (a) the energy per spin; (b) 
the order parameter Mq for MSW theory and for PEPS the peak height of the structure factor at 
Q = (it, it) (Neel) and Q = (q,q) (spiral); (c) the Gaussian spin stiffness; (d) the components of 
the spin-stiffness tensor (where p xx = p yy by symmetry, and pxy Ttial = 0). 
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The indication of a disordered phase drawn from the break down of MSW theory is 



further corroborated by the order parameter M [Fig. 15 (b)], which decreases strongly for 



J3/J1 0.39 - and for J3/J1 — > 0.52 + , and by the spin stiffness [Fig. 15 (c) and (d)], which 
is drastically reduced when approaching the above two limits. In particular, the Gaussian 
spin stiffness T is already strongly reduced for J3/J1 > 0.3. These results are consistent 
with the vanishing of the spin stiffness at J3/J1 = 0.35 that was found by ED of a system 
of 20 sites in Ref. [67]. 

A destabilization of magnetic order at around J3/J1 > 0.3 seems to be confirmed by the 



PEPS order parameter, Fig. 15 (b), which vanishes in the range 0.3 < J3/J1 < 0.5. Note 
that, again, we find that the PEPS order parameter deep in the Neel phase is similar to the 
MSW data, but that in the spiral phase MSW data for the order parameter lie well above 
the PEPS ones. 

In our calculations, despite using the same equations as in Ref. (30], we find a considerably 
larger breakdown region. However, the region where our calculations do not yield a result 
is very stable, i.e., it does not depend much on system size nor on the exact algorithm for 
solving the self-consistent MSW equations. 

The precise nature of the state in the candidate region for quantum-disordered behavior 
cannot be determined reliably by the use of MSW theory. From an analysis of the dimer- 
dimer correlations in the convergence regions, we can find no indications of any exotic 
disordered quantum state; on the contrary, PEPS results indicate a plaquette state in the 
region of maximal frustration J 3 w J\/2 



D. Ground state phase diagram of the J1J2J3 model 



1. MSW results 

After having investigated the two limiting cases of the J\ J2 and the J\ J3 models, we con- 
sider more generally the J1J2J3 model over the relevant parameter range < J2/J1, J3/J1 < 
1. As already seen in the case of the J x J 3 model, we observe a sizable parameter range over 
which the convergence of MSW theory breaks down, and which is then pointed out as a 
candidate region for non-magnetic behavior. We notice that, while convergence is achieved 
for any J2/J1 ratio at J 3 = 0, a region of convergence breakdown opens up by adding a 



26 




0.0 0.2 0.4 6 0.8 1.0 



J 2 IJ\ 



0.0 0.2 0.4 0.6 0.8 1.0 



c) 



1.0 
0.8 



^0.6 



Y 



0.4 
0.2 



0.0 




0.05 i.o 



•^partial 



0.0 0.2 4 0.6 0.8 1.0 








0.05 



0.0 0.2 0.4 0.6 0.8 1.0 



Figure 16: (a) Ground state energy per spin Eq, (b) order parameter Mq, (c) Gaussian spin 



stiffness T, and (d) Gaussian spin stiffness T partial calculated via Eq. (A10). Note that T and 
Impartial rise beyond the linear scale in the upper half of the plot. In the gray areas convergence of 
the self-consistent equations could not be reached. The blue lines are the classical phase boundaries. 

small J3 component around J2/J1 ~ 0.5. The energy per spin increases when approaching 



this region, showing the increased influence of frustration [Fig. 16 (a)]. The indications for 



a quantum disordered phase in the break-down region is corroborated by the decrease of the 



order parameter [Fig. 16 (b)] and the spin stiffness [Fig. 16 (c) and (d)] when approaching 



the break-down region. 

The nature of the phases where MSW reaches convergence can be seen in the ordering 



vector, which we display in Fig. 17 in comparison with the classical one. We find three 
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Figure 17: Ordering vector for the J1J2J3 model in linear color scale. In the gray area convergence 
of the self-consistent equations could not be reached, (a) x-component, and (b) y-component of the 
classical ordering vector; (c) x-component, and (d) y-component of the quantum mechanical MSW 
ordering vector. The blue lines are the classical phase boundaries. 

ordered phases: 1) For small J3/J1 and J2/J1 we find a Neel ordered phase. Its boundary is 
pushed upwards to higher values of J3/J1 with respect to the classical limit; 2) a columnar 
phase is found at small J3/J1 but larger J2/J1 > 0.6; 3) for large J3/J1 a spiral phase arises 
with an ordering vector Q = (q,q) that approaches Q = (7r/2,7r/2) for large J3/J1. As a 
consequence, a most dramatic effect of quantum fluctuations seems to be the disappearance 
of phase III in the classical phase diagram, characterized by magnetic order at a pitch vector 
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Figure 18: M (Q) [Eq.([5])] for a PEPS calculation on a 8 x 8 lattice with auxiliary dimension D = 3. 
A low value marks a destabilization of magnetic LRO. 

Qd = (?) 7f) with continuously varying g, in favor of the columnar phase and of a potentially 
quantum-disordered phase. 



2. Comparison to PEPS calculations 



In Fig. 18, we display the peak height of the static structure factor, Eq. ([5]), from a 
PEPS calculation on a 8 x 8 lattice with auxiliary dimension D = 3. We observe a broad 
asymmetric v-shaped region in which the magnetic order, quantified by the height of the 
peak in the structure factor, is strongly suppressed. We notice that this region is strongly 
reminiscent of (albeit broader than) the breakdown region of MSW theory. In particular, the 
asymmetry is due to the fact that the bottom of the "v" lies at J%j J\ > 0.5, a characteristic 
which is shared with the MSW phase diagram. While a thorough finite-size scaling analysis 
of the PEPS data would be necessary to determine the precise boundaries of the possible 
magnetically disordered regions, a quantitative information can be extracted even from the 
finite-size PEPS data concerning the location of the pitch vector of the dominant (long- 
ranged or short-ranged) magnetic correlations. 

Similarly to what happens in the above spin-wave calculations, a pronounced peak at the 
Neel ordering vector (7r, it) appears if both J 2 / ' J\ and J3/J1 are small, while at large J2/ J\ 
but small J3/J1 the structure factor is peaked at the columnar ordering vector (vr,0). For 
large J3/J1, finally, the peak is located at (q,q), where q tends to n/2. 
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Figure 19: Components of the ordering vector for a PEPS calculation on a 8 x 8 lattice with 
auxiliary dimension D = 3. 

IV. CONCLUSION 

In this work, we made use of Takahashi's modified spin-wave theory with ordering vector 
optimization to determine the ground state phase diagram of two paradigmatic models 
of two-dimensional frustrated antiferromagnetism: the S — 1/2 Heisenberg model on the 
SATL and on the J1J2J3 lattice. The optimization of the ordering vector shows dramatic 
quantum corrections to the ordering vector for spiraling states present in both models: such 
corrections show the general trend of promoting collinearly ordered states (either Neel or 
columnar states) against spiraling ones. Both for the triangular and the J1J2J3 lattice, 
MSW theory breaks down over a sizable region of parameter space, showing a dramatic 
suppression of the order parameter and of the spin stiffness as the breakdown region is 
approached: this finding is strongly suggestive of the appearance of quantum-disordered 
regions in the phase diagram of the models under investigation, an issue which is still under 
intense debate. The extent of the quantum-disordered regions estimated via MSW theory 
generally appears to be lower than that estimated by more accurate numerical techniques 
which take into account quantum fluctuations in a more complete fashion. Hence, one can 
draw two main conclusions from our results: on the one hand MSW might still converge 
to a magnetically ordered ground state even though the true ground state is disordered - 
although in this case it will probably feature a small value for the order parameter, or a 
small stiffness, suggesting that the magnetic order is not robust when dealing with quantum 
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fluctuations more accurately; on the other hand, the breakdown of MSW theory seems to 
be a strong indication that the true ground state is disordered. 

In particular, in the case of the SATL, MSW theory completely breaks down for suffi- 
ciently weak couplings between the chains composing the lattice, suggesting that the system 
remains in a disordered ID-like state even when the chains are coupled, as already predicted 
by recent variational approaches. A further disordered phase appears when the inter-chain 
couplings exceed the intra-chain ones: this phase is sandwiched in between the spiral phase 
of the nearly isotropic triangular lattice and the Neel phase appearing at large interchain 
couplings. In the case of the J1J2J3 lattice, a large breakdown region separates the Neel- 
ordered region for small J2 and J3, from the columnar-ordered region for J 2 > Ji/2 and small 
J3, and from the spiral phase at large J3. Hence, a general conclusion that we can draw 
from the study of these two models is that collinearly ordered phases (Neel and columnar) 
and spiral phases cannot be connected adiabatically - at least at the MSW level - but they 
are always separated by a breakdown region; this is a signal that in the true ground state 
collinear and spiral phases might always be divided by an intermediate quantum-disordered 
phase. 

Quantitative comparisons with more accurate methods (exact diagonalization, and vari- 
ational Ansatzes based on projected BCS states and projected entangled-pair states) reveal 
that MSW theory with ordering wave- vector optimization goes well beyond linear spin-wave 
theory in dealing with quantum effects, and it correctly accounts for the quantum correction 
to the ordering wave-vector of the ordered phases, and for the strong suppression (or total 
cancellation) of magnetic order in correspondence with the candidate regions for quantum- 
disordered behavior. Given its flexibility and its modest numerical cost, MSW theory serves 
therefore as a unique tool for the identification of novel quantum phases in strongly frustrated 
quantum Heisenberg antiferromagnets. 
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Appendix A: Modified spin-wave formalism for Heisenberg antiferromagnets 



In this Appendix, we shortly review the MSW formalism as applied to Heisenberg an- 
tiferromagnets. The full description of the approach - as applied to XY models - can be 
found in |31j . 

The Dyson-Maleev transformation |68| 169] maps the Heisenberg Hamiltonian, Eq. ([!]), 
to the non-linear bosonic Hamiltonian 



(hi) 



+ 



2S \ a\cij + diOjj — a\a}j(ijaj — a\aiaia^ (1 + cos (Q ■ r^)) 
2S \ a\cij + a^a^j — a^ajaj — a\aiaiaj (1 — cos (Q ■ r^)) 



(Al) 



+ 4 



S 2 — S (a\a,i + cijCij) + a\aia}-aj cos {Q ■ r^) + O 



S 



where <jj (a|) destroys (creates) a Dyson-Maleev boson at site i, S is the length of the spin, 
and Q the ordering vector. Here, we neglected the kinematic constraint which restricts 
the Dyson-Maleev-boson density n to the physical subspace n < 2S, given by the length 
of the spins S. Moreover, we dropped terms with six boson operators, which are of order 
0[n/ (2S) 3 ] and are negligible for n/(2S) < 1. Using Wick's theorem [70], and defining the 
correlators (a\aj) = F (ry) — |5y and (a>iCij) = = G( r ij)i the expectation value 

E = (7i) can be written as 



E 



2 ^ tij 

<i,i> 



S+--F(0) + F(r y .; 



S+--F(0) + G(Ty 



;i + cos(Q-r -)) 



[1 - C0S(Q • Tij)) 



(A2) 



After Fourier transforming, ay, = Yli a i e~ lk ' Vi , where N is the number of sites, and 
a subsequent Bogoliubov transformation, a& = cosh^a*. — sinh#fcat fe , and a_ k = 
— sinh 6*fc afc + cosh#/,a^ fc , we minimize the free energy under the constraint of vanishing 
magnetization at each site, (a|aj) = S [29]. (This guarantees that the kinematic constraint 
is satisfied in the mean.) This yields a set of self-consistent equations, 



tanh 29 k 



(A3) 
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with 

= a7 E ^ " cos & ' r v)) ^ (A4b) 
- (1 + cos (Q • rij )) Fij (1 - e**^)] - ii , 
where ji is the Lagrange multiplier for the constraint. The spin- wave spectrum reads 



^ = \JBl-Al. (A5) 

At T = 0, where n& = Vfc ^ 0, one finds that /i vanishes. This implies also the disap- 
pearance of the gap at k = that may exist for finite temperature. A vanishing gap is a 
necessary condition for magnetic LRO. It also enables Bose condensation in the k = mode. 
Separating out the contribution of the zero mode, {a\. =0 ak=o) /N = (a k=0 a k=Q ) /N = M 
(corresponding to the magnetic order parameter), one arrives at the zero-temperature equa- 
tions 

Fij = M + ± ^ cos (fc • r -) , (A6a) 

G y = M o + ^ E cos (fc ' r ^ ' (A6b) 

and the constraint of vanishing magnetization at each site becomes 

It is not a priori clear that the classical ordering vector Q cl correctly describes the LRO in 
the quantum system. To account for the competition between states with LRO at different 
ordering vectors Q we extend the MSW procedure by optimizing the free energy T with 
respect to the ordering vector Q. This yields two additional equations which must be added 
to the set of self-consistent equations, 

^=~X)^8in(QTy)r?.[4 + G?,.] =0, (A8a) 



dQ x 2 



(id) 



d 



F = ~\ E (Q • r y ) rj- [i* + G?-] = . (A8b) 



°Q, 2 (w> 
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In the SATL with NN interactions these simplify to Q y = and 



Q x = 2 arccos 



2F3 +GI 



(A9) 



where n = (1, 0) and r 2 = (1/2, ^3/2) are the lattice vectors. 



The values of and can now be calculated by solving self-consistently Eqs. (A4- A8). 
Through Wick's theorem the knowledge of the quantities F^ and allows the computation 
of the expectation value of any observable. 

a. Spin stiffness The optimization of the ordering vector allows a straightforward cal- 
culation of the spin stiffness, which gives a measure of how stiff magnetic LRO order is 
with respect to distortions of the ordering vector, and thus provides a fundamental self- 
consistency check of our approach. In fact, finding a small spin stiffness casts doubt on 
the reliability of the spin-wave approach in describing such a strongly fluctuating state, and 
hence suggests that the true ground state might be quantum disordered. 

The spin stiffness tensor is defined as p a p - 1 ' """" 



, ,., , evaluated at the optimized 

N dQadQfj q = qo' ^ 

ordering vector Q°. From this we can extract the parallel spin stiffness p\\ = | (p xx + p yy ) 
and the Gaussian spin stiffness T = det p. 

Since a change in Q affects the correlators F^ and G^, we must compute T self- 
consistently. After finding the optimal Q° by the self-consistent procedure described above, 
we calculate j^F (Q x , Qy) self-consistently for several fixed ordering vectors Q = Q° + AQ 
and fit a quadratic form to the results. Since the minimum in the free energy can be very 
shallow, this procedure can be affected by numerical noise. As an approximation to the true 
spin stiffness, the partial spin stiffness p^a tial can be computed via the partial derivatives, 
i.e., without recalculating the self-consistent equations. It reads 

^ - 'nScBq/ (A10) 



(id) 

We define T partial analogously to T as the determinant of the partial spin-stiffness tensor. 
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